Charge Inversion of Divalent Ionic Solutions in Silica Channels 
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Recent experiments (F.H.J. Van Der Heyden et al., PRL 96, 224502 (2006)) of streaming currents 
in silica nanochannels with divalent ions report charge inversion, i.e. interfacial charges attracting 
counterions in excess of their own nominal charge, in conflict with existing theoretical and simula- 
tion results. We reveal the mechanism of charge inversion by using all-atomic molecular dynamics 
simulations. Our results show excellent agreement with experiments, both qualitatively and quanti- 
tatively. We further discuss the implications of our study for the general problem of ionic correlations 
in solutions as well as in regards of the properties of silica-water interfaces. 

PACS numbers: 82.45.Mp,61.20.Qg,82.39.Wj 



I. INTRODUCTION 

Counterions in aqueous solution play a crucial role in 
the self-assembly of colloids and polymers, cell signaling, 
micronuidics and fuel cells. In recent years, there has 
been a considerable body of experimental and theoreti- 
cal work aimed at understanding the rich and complex 
phenomenology of ions in aqueous solution and their in- 
teractions with charged interfaces [J H, H, 0] • A repre- 
sentative example of this phenomenology is charge inver- 
sion (CI), where interfacial charges attract counterions 
in excess of their own nominal charge, thus leading to an 
interface whose effective surface charge is reversed. The 
first simulation study that showed evidence of CI with 
divalent ions was conducted by Torrie and Valleau [H , in 
which the authors used a primitive hard sphere model for 
the ions and a continuum solvent model with a dielectric 
constant. Still the origins of CI are poorly understood, as 
highlighted in a recent review article by Lyklema [6j. In 
particular, it is unclear to what extent existing theories 
explain the observed CI reported in experimental work. 

A recent experiment has provided a detailed investiga- 
tion of CI by multivalent counterions from an analysis of 
streaming currents, the electric currents resulting from 
applying a pressure gradient along a charged nanofluidic 
channel containing electrolyte solutions (see Fig. [T]) Q • 
The precision and accuracy of such experiments provide 
a perfect arena to test existing theoretical models of CI 
in particular and the more general problem of ion distri- 
butions. 

Different theories for CI have been proposed over the 
years. In Ref. § it was proposed that CI arises from the 
gain in correlation energy brought about by lateral cor- 
relations (LC) among mobile counterions. For divalent 
ions in silica nanochannels, this theory would predict CI 
at Cmv ~ 10 mM [3?|]-an order of magnitude lower than 
what is observed experimentally In Refs. H El it 
has been pointed out that the intrinsic discreteness of 
interfacial charges may give rise to CI by transverse cor- 
relations (TC), that is, correlations between counterions 
and interfacial charges. In particular, it has been shown 
that binding constants of divalent ions to oxygen atoms 
are typically of the order of Kb ~ 10 M _1 , thus imply- 
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FIG. 1: (Color online:) Cartoon picture of a silica nanochan- 
nel. Very near to the silica interface, we identify a first layer 
of bound counterions (layer 1), a layer of hydrated counteri- 
ons (layer 2) and a third layer (layer 3) , which corresponds to 
the diffuse layer and exhibits CI. 



ing CI at Ci nv ~ 1/Kb k 0.1 M. Numerical simulations 
with primitive models, where the solvent is considered 
implicitl y as a medium of dielectric constant e, have been 
reported [n], [12, [T^j , and show general agreement with LC 
theories, thus implying lower than the Cj n „ s=s 400 mM 
observed in experiments Q- Other theoretical work has 
focused in the computation of mobilities [l4[ as a function 
of (^-potential, but as reviewed in [6j], the interpretation 
of experimentally determined (^-potentials is not entirely 
clear. A clear understanding of the CI phenomenon may 
therefore require approaches where water is considered 
explicitly and the interface is modelled realistically, with 
more structure than a uniform surface charge. As dis- 



cussed in [10|, |l5|, these effects are expected to play a 
critical role. 

In this paper, we investigate the origin of CI with 
atomic resolution by using molecular dynamics (MD) 
simulations of a charged silica interface in contact with 
an aqueous solution with divalent ions. This study 
is further motivated by recent interest in understand- 
ing the structure of pure water in contact with silica 
flo . IT7L [l8l Il9l . [20j , as some evidence suggests the pres- 
ence of highly ordered interfacial water structures. It has 
not been investigated, however, how divalent ions, which 
have a very cohesive hydration sheath, may modify inter- 
facial water. In addition, our results have implications in 
the field of microfiuidics because silica is the most com- 



mon material used in nano- and microfluidic devices 2l| , 
including the recently synthesized nanoporous function- 
alized silica thin films [22|, l23| that are being investigated 
as desalination membranes 12411. 



II. SIMULATION DETAILS 

We consider the system schematically shown in Fig. [1] 
where an interface of amorphous silica is in contact with 
an aqueous solution. Throughout this letter, we will con- 
sider ionic strengths such that Xd << h, where Ad is the 
Debye length and h is the separation between the two sil- 
ica interfaces (see FigfT]). We can therefore neglect the 
interactions between the two silica interfaces and perform 
our MD simulations on a single interface. 

Amorphous silica was generated according to a previ- 
ously reported procedure [III [26[ . We start with a bulk 
a-quartz crystalline substrate that is heated to high tem- 
perature (~ 2500 K) and quenched to 300 K. We then 
generate two free surfaces in the z-dimension and further 
anneal the substrate until the density of bond defects 
on each interface is representative of the experimen tally 
observed density on amorphous silica (4.0-5.0/nm 2 )[27|]. 
The silicon and oxygen defects are terminated with - 
OH and -H respectively. The silica interface is given a 
negative charge by randomly selecting -OH groups and 
removing the corresponding proton, leaving a net sur- 
face charge of ~ le/110 A~ 2 , in line with experimental 
results 7]. The amorphous silica had a total area of 
68.7076 x 68.082 A 2 in contact with an aqueous solution 
containing water (~ 8500 water molecules) and divalent 
ions. Water was modelled explicitly by using the TIP3P 
model [28] , the two cations used in the simulation, (Ca 2+ , 
Mg 2+ ), were modeled with Aqvist's parameters [29| and 
Cl~ with the OPLS forcefield [3(J. The silica was mod- 
eled by the OPLS force field described in (3l| . 

The simulations reported in this article were performed 
with the LAMMPS molecular dynamics code [32j. In all 
simulations, the temperature is fixed at 300 K and 
controlled with a Nose-Hoover thermostat and a time 
constant of 0.01 fe _1 . The SHAKE algorithm [33[ was 
used to constrain the bond lengths and angles of the wa- 
ter molecule and the bond lengths of the OH bonds on 
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FIG. 2: (Color online:) Evolution of the Ca 2+ charge (blue 
circle) and CI - charge (red square) relative to the interfacial 
charge as a function of time for the different layers. Layer 1 
represents bound oxygens (defined such that 0-Ca 2+ distance 
is less than 3 A) layer 2 is defined by the size of the hydrated 
divalent ions (0-Ca 2+ distance is within 3 and 6 A) and 
layer 3 is the last layer before bulk values are attained, with 
0-Ca 2+ distance < 12 A. (see Fig.[T]for a schematic definition 
of the different layers). Results are for 1.0 M concentrations. 



TABLE I: Relative charge population of the ionic species of 
Ca 2+ and Cl~ in each layer. Values in parentheses represent 
fluctuations in the last digit (s). 



System 


1 st layer 
Ca 2+ CP 


2 nd layer 
Ca 2+ cr 


3 rd layer 
Ca 2+ cr 


0.5 M 
1.0 M 
1.5 M 


0.68(41) 0.02(1) 
0.99(29) 0.10(2) 
0.92(40) 0.11(1) 


0.54(23) 0.13(2) 
0.69(13) 0.43(6) 
0.85(27) 0.55(19) 


0.33(7) 0.36(8) 
0.70(18) 0.88(12) 
1.03(26) 1.14(17) 



the silica substrate surface. A 1 fs time step was used and 
the van der Waals (vdW) interaction is cut off at 10 A. 
A slab version of the particle-particle particle-mesh al- 
gorithm [U is used to compute long-range Coulombic 
interactions. An average run lasted around 12 ns, with 
data being collected and analyzed over the last 8 ns. A 
typical run took 4608 node-hours per 12 ns simulation. 



III. SIMULATION RESULTS 

It is convenient to divide the region near the silica in- 
terface into three layers as depicted in Fig. [TJ The first 
layer is defined so that it includes all partially dehydrated 
divalent ions bound to silica oxygens. The second layer 
includes all counterions within a distance of 6 A (mea- 
sured from the center of the silica oxygens) and basi- 
cally consists of hydrated counterions. The third layer 
extends up to a region of 12 A, beyond which bulk den- 
sity profiles are attained. We recall that in textbooks in 
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FIG. 3: (Color online:) Net charge density on the third layer 
(see Fig. [1] for a definition of the layers) for different concen- 
trations of CaCb- The inset shows the profile of the charge 
density. 
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FIG. 4: (Color online:) Radials distribution function of depro- 
tonated silica oxygen atoms (Obare) and protonated oxygen 
atoms (OH) with Ca 2+ . The inset shows the rdf for (Ob are )- 
divalent ions at a larger vertical scale. 



physical chemistry [35| , layer 1 is referred to as the Inner 
Helmholtz plane, layer 2 as the Outer Helmholtz plane 
and layer three as the diffuse layer. The time evolution 
for the population of both Ca 2+ and Cl~ on the different 
layers are shown in Fig. [2] at 1.0 M concentration. The 
population is expressed as a function of total charge of 
the ionic species within the layer relative to the total net 
interfacial charge. As is clear from Fig. ^ the popula- 
tion of the different layers reaches equilibrium after an 
initial period of about 3 ns. We further measured res- 
idence times and found that the average residence time 
for Ca 2+ ions were ~ 2.0 ns, ~ 200 ps and ~ 30 ps in 
the first, second and third layers respectively. The to- 
tal simulation time for our simulations was 12 ns thus 
providing ample time to compute equilibrium quantities. 
No significant differences were observed for other CaCl2 
concentrations. A summary of the average populations 
within different layers is shown in Table [IJ 

Results for MgCl2, however, show significant differ- 
ences, as residence times in the first and second layer 
are surprisingly long. As an example, within 12 ns, only 
a single Mg 2+ ion was exchanged from the first to the 
second layer. It is therefore possible that the first and 
second layer for MgCl2 are not thermalized. The results 
for the MgCl2 solutions are also of interest as will be 
clear, so in view of these caveats we will relegate a de- 
tailed summary of these results to Appendix [Al 

The analysis of the time series Fig. [2] and Table U shows 
that the third layer has a net negative charge, that is, 
there is an excess of Cl~ charge over divalent charge, 
which is the signature of CI. This is clear from Fig. [31 
where the charge density as a function of distance from 
the interface is shown. The results exhibit no CI at 0.2 M 
(not shown) and marginal CI for 0.5 M CaCl2 , with CI 
considerably amplified at larger concentrations. These 



results show conclusive evidence for a Stern layer con- 
sisting of divalent ions (first and second layers) followed 
by an "inverted" (negatively charged) diffuse layer. In 
all cases, charge densities reach bulk values (within the 
accuracy of our simulations) at 12 A from the interface. 

The structure of the ions near the silica interface is an- 
alyzed using radial distribution functions (rdfs), shown 
in Fig. H] for 1.0 M solutions. The rdf corresponding to 
divalent ions with deprotonated silica oxygens (Ohare) 
show sharp peaks at a distance roughly equivalent to the 
sum of the crystallographic radius of oxygen and diva- 
lent ions providing clear evidence for binding of partially 
dehydrated counterions to oxygens. The inset of Fig. [4] 
shows a much stronger binding for Ca 2+ ions, sufficient to 
virtually neutralize the interfacial charge, as clear from 
Table [I] The rdf also shows a maxima roughly corre- 
sponding to the sum of the radii of the hydrated divalent 
ion and an oxygen. Although not bound, hydrated ions 
are strongly correlated with silica oxygens, an effect that 
becomes more dramatic with Mg 2+ ions. 

We further investigated the structure and orientation 
of water near the interface and its relation to the observed 
ion distributions. A detailed study of pure water near 
amorphous silica interfaces will be reported elsewhere 
[3^ | . here we just briefly discuss the most salient aspects 
related to ion distributions. It is found that the average 
lateral separation of water molecules near charged silica 
are identical with pure water results, but average dis- 
tances between the water molecules and oxygens on the 
silica interface are significantly altered (4.0 A as opposed 
to 3.2 A for neutral silica), which we interpret as be- 
ing related to the Ca 2+ hydration sheaths pushing away 
other water molecules. Our MD simulations do not show 
evidence for any effect on the counterions attributable to 
a pre-existing interfacial water ordering. 
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We briefly discuss the properties of the bulk solution as 
the concentrations discussed in this paper cannot be con- 
sidered as dilute, and activity coefficients deviate consid- 
erably from unity, ft is found that the bulk rdf of Ca 2+ - 
Cl _ ions shows strong correlations with a significant peak 
at a distance corresponding to the sum of the crystal- 
lographic radius of both ions (results not shown), thus 
providing clear evidence for Bjerrum pairing. A quan- 
titative analysis shows that the bulk solutions contain 
8(2)%, 17(5)%, and 25(7)% of Bjerrum pairs for concen- 
trations of 0.5 M, 1.0 M and 1.5 M respectively. 

IV. DISCUSSION AND CONCLUSIONS 

We now estimate the magnitude of the streaming cur- 
rent that would be measured in the systems investigated. 
We consider a channel with length L, width w and height 
h (see Fig. [T]). The streaming current is 

l-h/2 

Istr — W dzp c {z)v{z) (1) 

J-h/2 

where p c (z) is the net charge at position z and v(z) is 
the velocity of the fluid along the channel. This velocity 
is given by Poiseuille flow v(z) — ^^{z 2 — (f) 2 ), where 
r\ ~ 0.9 cP is the shear viscosity of water. Using the 
values of p c {z) from our simulation, explicit values for the 
streaming current are obtained. For wide silica channels 
(h » Xd), this expression is further related to the £- 
potential 

„ APh f°° , , . . whe r , 

Istr = 2W—— / dzzp c { Z ) = -AP- U 2 

2r]L J AirrjL 

where p c {z) is the net charge distribution on an infinite 
slab and e r ~ 80 is the dielectric constant of water. The 
value of the ^-potential is defined at the plane where the 
fluid velocity is zero, which we assume coincides with 
the boundary between the Stern and the diffuse layer Q . 
Using h — 140 nm,w =50 fim and a length of the channel 
L = 4.5 mm as typical values [7J, we obtain I str ~ 1 
p A/bar (£ w 30 mV) at concentrations of 1 and 1.5 M 
(at 0.5 M, < I str < 0.3 pA/bar). 

This paper has presented all-atomic MD simulations 
of divalent ionic solutions in silica channels. Our simu- 
lations show no CI at 0.2 M so 0.2 < c inv < 0.5 M. At 
concentrations larger than 0.5 M, CI is substantial, both 
for CaCl2 and MgCl2. At concentrations of 1.0 M, we 
calculate a streaming current of the order of 1 pA/bar. 
These results compare extremely well with the experi- 
ments in Ref. [7|. 

Our simulations provide a detailed understanding of 
the structure of the ions near the interface and the under- 
lying mechanism for CI. It is shown that the Stern layer 
consists of both bound, partially dehydrated (layer 1 or 
Inner Helmholtz plane) , and mobile hydrated counterions 
(layer 2 or Outer Helmholtz plane). Counterions in layer 
2 are strongly localized near the silica oxygens. In CaCl2 



TABLE II: Relative charge population of the ionic species of 
Mg 2+ and Cl~ in each layer. Values in parentheses represent 
fluctuations in the last digit (s). 

System 1 1 st layer 1 2 nd layer 1 3 rd layer 

M g 2+ cr M g 2+ cr M g 2+ cr 

0.5 M 0.26 0.03(1) 0.99(20) 0.14(2) 0.32(6) 0.38(8) 
1.0 M 0.15 0.08(1) 1.54(39) 0.43(11) 0.65(8) 0.84(13) 



solutions, the population of the Inner Helmholtz Plane 
is significantly larger than the Outer Plane. Overall, our 
results support a picture where CI is due to electrostatic 
correlations between counterions and discrete interfacial 
groups (TC) or Bjerrum pairing [9j. This is particularly 
clear in ions with softer hydration sheaths, such as Ca 2+ . 

In summary, in this paper we provided a detailed ac- 
count of the mechanism of CI of divalent aqueous so- 
lutions near silica, in excellent agreement with experi- 
mental results Q. By conducting simulations in which 
the water is explicitly modeled and the silica substrate is 
modeled in a realistic manner with discrete charges, we 
are able to capture the atomistic physics (i.e. structure 
of the ions and water near the charged species on the 
silica substrate) that plays such a large role in interfacial 
phenomena like that present in this problem. It remains 
as a future challenge to develop simpler analytical models 
and further quantitative experimental tests. 
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APPENDIX A 

Table |TT] shows the average populations of the Mg 2+ 
and Cl~ ions in the three layers as defined in Fig.[TJ The 
values in Table |TT] show that the first two layers have a 
net positive charge, and the third layer has a net negative 
charge. This behavior is a signature of charge inversion 
and it is the same general behavior that was observed 
in the CaCl2 systems. The difference between the two 
systems is that in the MgCl2 systems the silica charge 
is compensated by Mg 2+ in the second layer, whereas in 
the CaCl2 systems the silica charge is compensated by 
Ca 2+ in the first layer. 

Charge inversion in the MgCl2 systems is also obvious 
in Fig. [5l where the charge density as a function of dis- 
tance from the interface is shown. The data from the 
CaCl2 systems are also plotted in Fig. [5] so that they can 
be compared to the MgCi2 systems. In both cases, there 
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FIG. 5: (Color online:) Net charge density on the third layer 
(see Fig. [1] for a definition of the layers) for different concen- 
trations of MgCk and CaCL> ■ The inset shows the profile of 
the charge density. 



is marginal CI of nearly identical magnitude at a concen- 
tration of 0.5 M. Then at a concentration of 1.0 M, both 
systems show significantly larger CI. These results for the 
MgCl2 systems show conclusive evidence for a Stern layer 
consisting of Mg 2+ ions (first and second layers) followed 
by an "inverted" (negatively charged) diffuse layer. 

In MgCl2 solutions, the number of divalent ions within 
the Stern layer is similar as in CaCl2 but its structure is 
quite different, as there are very few bound counterions. 
We raised the possibility that this may reflect that the 
Stern layer may not be thermalized, and that longer sim- 
ulations (at least ~ 1 fis) might reveal a structure more 
similar to the CaCb case. We point out, however, that in 
Ref. [§] it was concluded that binding constants of Mg 2+ 
to oxygen atoms are generally small, a result attributed 
to the compactness of the Mg 2+ hydration sheath. In any 
case, whether thermalized or not, our simulations show 
CI and, as apparent from the rdf in Fig. 01 correlations 
between Mg 2+ and charged interfacial groups are shown 
to play a fundamental role. 
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